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1. Introduction 

A monochromatic wave u propagating in a heterogeneous medium is governed by the Helmholtz 
equation 

(1) AM(r) + a;^(l + u(r))M(r) = 0, r G d=2,3 

where u G C describes the medium heterogeneities. For simplicity, we choose the physical units 
such that the wave velocity is unity and the wavenumber equals the frequency cu. 

The data used for imaging is the scattered field = u — v} governed by 

(2) Art® + = —oj^yu 


or equivalently the Lippmann-Schwinger integral equation: 

(3) n®(r) = oj'^ [ u(r') (n‘(r') + u®(r'))G(r,r')dr'. 

Jr3 


Here 


( 4 ) 


G(r,r') 


gia;|r —r'l 

47r|r—r^l ’ 


iH, 


d = 3 
d = 2 


is the Green function for the background propagator (A + uP‘)~^ where is the zeroth order 
Hankel function of the first kind. 

We consider two far-held imaging geometries: paraxial and scattering. In the former, both the 
object plane and the image plane are orthogonal to the optical axis while in the latter emission and 
detection of light can take any directions. In the former, we take u® as the measured data and in 
the latter we take the scattering amplitudes (see 0 below) as the measured data. 

• Paraxial geometry: For simplicity, let us state the 2D version. Let {z = zq} be the object 
line and {z = 0} the image line. With r = (x, ZQ),r' = [x', 0), we have 


( 5 ) 


\x,zq) = / u ( x ', 0 ) + 
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Figure 1. Two imaging geometries: (a) Diffraction (b) Scattering. 


where C is a complex number. 

• Scattering geometry: The scattered field has the far-held asymptotic (Born and Wolf 
1999) 

( - / 1 W r 

(6) = + f-=^, d = 2,3 

where the scattering amplitude A has the dimension-independent form 

(7) vl(r,d) = ‘^ [ u{r')(u\r') + u^{r'))e-^‘^^'-^dr'. 

jRd 


Note that since tt in ([^ and Q is part of the unknown due to multiple scattering, the inverse 
problem is a nonlinear one. To deal with multiple scattering effects in compressive sensing, it is 
natural to split the inverse problem into two stages: In the hrst stage we recover the masked objects 

V{x) = i/(x, 0) 0)-b 0)) (paraxial geometry) 

F(r) = i/(r) (n'(r) -b u®(r)) , (scattering geometry) 

with the Fourier-like integrals in Q and 0 as the sensing operators. In the second stage, we 
recover the true objects from the masked objects. 












For the most part of the article, however, we will focus on the first stage or make the Born 
approximation to linearize the imaging problem and turn to the multiple scattering effect only in 
Section [9l 


2. Outline 


In Section]^ we review the basic elements of compressive sensing theory including basis pursuit 
and greedy algorithms (orthogonal matching pursuit, in particular). We place greater emphasis on 
the incoherence properties than on the restricted isometry property because the former is much 
easier to estimate than the latter, even though the latter can also be established in several settings 
as we will see throughout this article. One thing to keep in mind about incoherence is that it is far 
beyond the standard notion of coherence parameter, which is the worst case metric (see © below). 
The incoherence properties are fully expressed in the Gram matrix of the sensing matrix, also 
known as the coherence pattern. Second thing noteworthy about incoherence is that the standard 
performance guarantees expressed in terms of the coherence parameter often underestimate the 
actual performance of algorithms. Its usefulness primarily lies in providing a guideline for designing 
measurement schemes. 

In Section we consider the Fresnel diffraction with the pixel basis. The pixel basis, having a 
finite, definite size, is emphatically not suitable for point-like objects. Indeed, in order to build 
incoherence in the sensing matrix, it is imperative that the wavelength be shorter than the grid 
spacing. In other words, the pixel basis is suitable only for objects that are decomposable into 
“smooth” parts relative to the wavelength. The sparsity priors then come in two kinds: (i) there 
are few such parts with 1-norm as proxy (ii) there are few changes from part to part with the total 


variation as proxy (Section 4.1). In the context of Fourier measurement, we introduce the notion 


of constrained joint sparsity to connect these two sparse priors and discuss basis pursuit (Section 

4.2) and orthogonal matching pursuit for joint sparsity (Section |4.3[ ). 

In contrast to the pixelated objects, point objects naturally do not live on grids. Such a problem 

arises in applications e.g. discrete spectral estimation among others. There is this fundamental 

tradeoff in using a grid to image point objects with the standard theory of compressive sensing: 

the finer the grid, the better the point objects are captured but the worse the coherence parameter 

becomes. In Section we use the notion of coherence band to analyze the coherence pattern and 

design new compressive sensing algorithms for imaging well separated, off-grid point objects. In 
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addition to off-grid point objects, the coherence-band techniques are also useful for imaging objects 
that admit a sparse representation in highly redundant dictionaries. One celebrated example is the 


single-pixel camera discussed briefly in Section 5.4 


In Section we discuss Fresnel diffraction with sparse representation in the Littlewood-Paley 
basis which is a slowly decaying wavelet basis in stark contrast to the pixel basis and the point-like 
objects. In this basis, the sensing matrix has a hierarchical structures completely decoupled over 
different scales. In Section we discuss near-held diffraction in terms of angular spectrum which 
works out nicely with the Fourier basis. 

In Section]^ we consider inverse scattering with the pixelated as well as point objects. Here we 


focus on the design of sampling schemes (Section 8.2) and various coherence bounds for different 
schemes (Section |8.3| ). 

In Section we discuss multiple scattering of point objects and the appropriate techniques for 
solving the nonlinear inverse problem. The keys are the combination of the coherence-band and 
the joint sparsity techniques developed earlier. 

In Section we discuss inverse scattering with extended objects sparsely represented in the 


Zernike basis. In Section 11 we discuss interferometry with incoherent sources in astronomy. As 
a consequence of the celebrated Van Citter-Zernike theorem, the resulting sensing matrix has a 
similar structure to that for scattering with multiple inputs and outputs. The difference between 
them lies in the fact that for interferometry the inputs and outputs are necessarily correlated while 
for scattering the inputs and outputs can be independent. As a result, the (in)coherence properties 
of interferometry are more subtle and it is an ongoing problem to search for the optimal sensor 
arrays in optical interferometry in astronomy. 


3. Review of compressive sensing 

A distinctive advantage of compressive sensing is accounting for the finite, discrete nature of 

measurement by appropriately discretizing the object domain. 
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By a slight abuse of notation, we use || • ||p to denote the p-norm {p > 1) of functions as well as 
vectors, i.e. 

(8) ||/||p = Q IfirWdrY", f€LPiR<^) 

\ i/p 

l/.-rj , 

and ||f||o (the sparsity) denotes the number of nonzero components in a vector f. 

By discretizing the right hand side of ([^ or Q and selecting a discrete set of data on the left 
hand side, we shall rewrite the continuous models in the form of linear inversion 


(9) 


IfIL = 


N 

E 


( 10 ) 


g = + e 


where the error vector e G is the sum of the external noise n and the discretization error d due 
to model mismatch. By dehnition, the discretization error d is given by 


( 11 ) 


d = g — n — . 


Consider the principle of basis pursuit denoising (BPDN) 


( 12 ) 


min ||h||i, s.t. ||g — ^h ||2 < ||e ||2 = e. 


When e = 0, (12) is called basis pursuit (BP). With the right choice of the parameter A, BPDN is 


equivalent to the unconstrained convex program called the Lasso (Tibshirani 1996) 


(13) 


min -||g - $z ||2 + Ae||z||i. 

z 2 


Both BPDN (12) and Lasso (13) are convex programs and have numerically efficient solvers (Chen 
et al. 2001, Boyd and Vandenberghe 2004, Brucskstein et al. 2009). 

A fundamental notion in compressed sensing under which BP yields a unique exact solution is 
the restrictive isometry property (RIP) due to Candes and Tao 2005. Precisely, let the restricted 
isometry constant (RIC) 6s be the smallest nonnegative number such that the inequality 


k {1 - (5s)||h||2 < ll^hlli < k( 1 + 5s)||h||2 
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holds for all h G C'^ of sparsity at most s and some constant k > 0. RIP means a sufficiently small 
62 s (see (14) below). 


Now we recall a standard performance guarantee under RIP. 


Theorem 1. (Candes 2008) Suppose the RIC of ^ satisfies the inequality 
(14) 82s<V2-1 


with K = l. Then the solution R of BPDN (12) satisfies 


(15) 


IR-flb < Cis-^/^Wi - &\\i + C2e 


for some constants Ci and C 2 where consists of the s largest components, in magnitude, off. 


Remark 1. For general n 1, we consider the normalized version of (10) 



and obtain from (15) that 


(16) 


R-f||2 < 

VK 


Note however that neither BPDN or Lasso is an algorithm by itself and there are many different 
algorithms for solving these convex programs. Some solvers are available on-line, e.g. YALLl and 
the open source code Ll-MAGIC (http://users.ece.gatech.edu/~justin/llmagic/). 

Besides convex programs, greedy algorithms are an alternative approach to sparse recovery. A 
widely known greedy algorithm is the Orthogonal Matching Pursuit (OMP) (Davis et al. 1997, 
Pati et al. 1993). 
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Algorithm 1. Orthogonal Matching Pursuit (OMP) 

Input: ^,g. 

Initialization: = 0, = g and = 0 

Iteration: For j = 1,s 

1) ima.^ = argmaxj | \ ,i ^ 

2) = 5^-1 U { i ^^} 

3) P = argminh ||$h — g ||2 s.t. supp(h) CS^ 

4) r-^ = g — 

Output: P. 


OMP has a performance guarantee in terms of the coherence parameter dehned by 

(17) = jj^f^ 

where is the A:-th column of is the pairwise coherence parameter and the totality 

[fi{k,l)] is the coherence pattern of the sensing matrix Here and below f denotes the conjugate 
transpose. 


Theorem 2. (Donoho et al. 2006) Suppose that the sparsity s of the signal vector f satisfies 
(18) /r($)(2s-l) + 2& < 1 

/min 


where /min = min \fk\- Denote by f*, the output of the OMP reconstruction. Then 
k 

(a) f* has the correct support, i.e. supp{p) = supp{P where supp{{) is the support o/f. 

(b) f* approximates the object vector in the sense that 


(19) 


|f* - fib < 


pi + p,- IJ,S 


Incoherence or RIP often requires randomness in the sensing matrix which can come from the 


randomness in sampling as well as in illumination. Between the two metrics, incoherence is far 
more flexible and easier to verify for a given sensing matrix. However, performance guarantees in 
terms of the coherence parameter such as (18) of Theorem tend to be conservative. 










4. Fresnel diffraction with pixel basis 


As a first example, we consider the imaging equation for Fresnel diffraction. We shall write 
Q in the discrete form (10) by discretizing the right hand side of ([^ and selecting a discrete set 
of scattered field data for the left hand side. 

We approximate the masked object 


f20i 


by the discrete sum on the scale i 

N 

(21) Vi{x) = ^Kj- k)Viek), V{ik) = i^{£k)u{ek, 

k=l 

where 


( 22 ) 


b{x) = 


1 , x£[-k,k] 

0, else. 


is the localized pixel “basis”. We assume that is a good approximation of the masked object for 
sufficiently small i in the sense lim£_,.o ||F — V^||i = 0. 

Moreover, we assume that is sparse in the sense that relatively few components V{ki) are 
significant compared to the number of grid points N. Note that sparse objects in the pixel basis are 
not point-like. Point objects typically induce large gridding errors and requires techniques beyond 
standard compressive sensing reviewed in Section (cf. Section . 

To proceed, we shall make the Born approximation and set u^{x,0) = 1 (i.e. normal incidence 
of plane wave). 

Let Xj,j = be the sampling points on the image/sensor line and define 

. 

Set the discretized, unknown vector f G as 


k = l, ...,N 
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and the data vector g G as 


^ U%X,,Zo) ^.i^^2/^2zo) . ^ ^ ^ 

^ Cibi^j) 


where 

(24) 


6(0 = / 6(x)e-*2-^«dx = 


sin (ttO 

TT^ 


As a result, @ can be expressed as (|10[) with the sensing matrix 


(25) 


$ = 


$1 . . . 4>7v 


E C^xAf^ ^ 


^—2‘Ki^jk 


M 


k = l,...,iV. 


J j=i 


A sensing matrix whose columns have the same 2-norm (as in (25)) tends to enjoy better perfor¬ 
mance in compressive sensing reconstruction. 


When 0 are independent uniform random variables on [—1/2,1/2], (25) is the celebrated random 
partial Fourier matrix which is among a few examples with a relatively sharp bound on the RIP 
given below. 


Theorem 3. (Rauhut 2008) Suppose 

M 


(26) 


1 


, ^ khr k\n N \n -, e E (0,1) 

InM “ e ^ ^ 


for given sparsity k where c is an absolute constant. Then the restricted isometry constant of the 

6 k <6 


matrix (25) satisfies the bound 


with probability at least 1 — e. 


Remark 2. To apply Theorem^ in the context of Theorem\^ we can set k = 2s and 6 = \/2 — 1. 


Ineq. (26) then implies that it would take roughly 0{s), modulo some logarithmic factors, amount 


of measurement data for BPDN to succeed in the sense of (15). 

On the other hand, the coherence parameter p, typically scales as as we will see in 

Theorem^ so, in view of the condition (18) in Theorem^ the amount of needed data is 0(s^}, 
significantly larger than 0{s) for 1 <C s <C A^. 
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While this observation is usually valid in the case of OMP, it needs not apply to other greedy 
algorithms such as Subspace Pursuit (BP) whose performance guarantee requires 0{s), up to loga¬ 
rithmic factor, amount of data (Dai and Milenkovic 2009). 


The fact that are independent uniform random variables on [—1/2,1/2] implies that Xj are 
independent uniform random variables on [—A/2, A/2] with 

27rzo 


(27) 


A = 


ojI 


in view of (23). Viewing ^ as the resolution length of the imaging set-up we obtain the resolution 
criterion 


(28) 


£ = 


27rzo 

Auj 


which is equivalent to the classical Abbe or Rayleigh criterion. 

Now let us estimate the discretization error vector d in ( |11[ ). Define the transformation 7” by 

1 


(TV),- = 






cf. 0- By definition 


d = TV- TVe 


we have 
(29) 


|tl||oO ^ 


IR- 


m = 


sin (tt^) 


imiuj \b{f,j)\' 

For ^ G [—1/2,1/2], min |S(.^)| = 2/7r and max |6(^)| = 1. Hence 

iry/M 


dh < ||d||oo\/M< 


|d||2 ^ 7rCVM\\V-Vi\\i 


(30) 

and 


which can be made arbitrarily small by setting i sufficiently small while holding M fixed and 


maintaining the relation (28). 
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Figure 2. The original 256 x 256 Shepp-Logan phantom (left), the Shepp-Logan 
phantom and the magnitudes of its gradient with sparsity s = 2184 (Fannjiang 2013. 
Reprinted with permission). 


4.1. Total variation minimization. If the masked object V is better approximated by a piecewise 
(beyond the scale i) constant function then the sparsity prior can be enforced by the discrete 
total variation 


|h||tv = ^|A/i(j)|, Ah{j) = hj+i - hj. 


Instead of (12) we consider a different convex program, called total variation minimization (TV- 
min) 


(31) 


min ||h||tv, s.t. ||g — ^h ||2 < e. 


cf. (Candes et al. 2006, Rudin and Osher 1994, Rudin et al. 1992, Chambolle 2004, Chambolle 
and Lions 1997). 

For two-dimensional objects h{i,j),i,j = l,...,n, let h = (hp) be the vectorized version with 
index p = j + {i — l)n. The 2D discrete (isotropic) total variation is given by 

l|h||tv = V|Ai/i(z,j)|2 -h |A2/i(i,j)P, 

Ai/i(f,j) = (/^(* + l>j) - ^(bi), ^2Hi,j) = h{i,j + l)-h{i,j)) . 
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Figure 3. BPDN reconstruction without external noise (left) and TV-min recon¬ 
struction with 5% noise (right) (Fannjiang 2013. Reprinted with permission). 


Fig. I and Fig. are a numerical demonstration of TV-min reconstruction of 2D object (the 
phantom). Fig. [^shows the original image and its gradient which is sparse compared to the original 
dimensionality. Fig. shows the reconstruction with BPDN (left) and TV-min (right). TV-min 
performs well as expected because the TV-sparsity is the correct prior for the object. On the other 
hand, BPDN performs poorly because the Ll-sparsity is the wrong prior. 


4.2. BPDN for joint sparsity. The close relationship between (31) and (12) can be seen from 
the following equation for the ID setting 


(e2-«^-l)5, = - fk). 

k 

In other words, the new data vector g = ((e^’^*^^ — l)gj), the new noise vector e = ((e^’^*^-’' — l)ej) 
and the new object vector f = {fk+i — fk) are related via the same sensing matrix as for BPDN. 
Clearly, \ej\ < 2\ej\,j = 1,...,M. Moreover, if ej are independently and identically distributed, 
then ej are also independently and identically distributed with variance 

E|ejp = - 1^ x Eje^f = 2E\ej\‘^ 


when is the uniform random variable over [—1/2,1/2]. Hence for large M the new noise magnitude 

||e ||2 ~ -v/2||e||2. Here and below E denotes the expected value 
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The similar relationship exists in the 2D case. Let fj = Ajf which satisfy the linear constraint 
(32) Alfa = Aafi. 

Dehne 


gi = - l)gj], g 2 = - l)gj] 

ei = - l)ej], ea = - l)ej] 


where ^j,gj,j = are independent uniform random variables over [—1/2,1/2]. Then F = 

[fi,fa] G G = [gi,g 2 ] G and E = [ei,ea] are related through 

G = [$fi,$f2] +E 


subject to the linear constraint (32). This formulation calls for the L^-minimization (Fannjiang 
2013) 


(33) 


min ||[hi,ha]|| 2 ,i, s.t. ||G - [$hi, $ha]||F < ||E||f, 


subject to the constraint 


(34) 


Aahi = Aiha 


where || • ||f is the Frobenius norm and || • Hap is the the mixed (2, l)-norm (Benedek and Panzone 
1961, Kowalski 2009). 


(35) 


|X||a,i = ^||rowj(X) 


The reason for minimizing the mixed (2, l)-norm in (33) is that fi and fa share the same sparsity 
pattern which should be enforced. 

To get a more clear idea about ||E||f, we apply the same analysis as above and obtain 


|ei||i ~ E||ei||2 = 2E||e||2, i = l,2. 


for sufficiently large M. 
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The convex program (33)-(34) is an example of BPDN with constrained joint sparsity. More 


generally, suppose that the columns of the unknown multi-vectors F G share the same support 

and are related to the data multi-vectors G G noise multi-vectors E G via 

(36) G = $2f2, ^jfj] + E 

subject to the linear constraint £F = 0. 

For this setting, the following formulation of BPDN with joint sparsity is natural 

(37) min||H|| 2 ,i, s.t. ||G - [$ihi, $ 2 h 2 ,..., $,/hj] ||f < e, s.t £H = 0, 
with e = ||E||f. 

4.3. OMP for joint sparsity. Next we present an algorithmic extension of OMP for joint-sparsity 
(Cotter et al. 2005, Chen and Hua 2006, Tropp et al. 2006)to the setting with multiple sensing 
matrices (|36|) (Fannjiang 2013). 


Algorithm 2. OMP for joint sparsity 

Input: {$j},g,e > 0 

Initialization: = 0, = G and ^ 

Iteration: For k = 1, 2, 3, • • • 

1) *max = argmaxj X]/=i where is the conjugate transpose of i-th column of 

2) 5^ = 5^-1 U 

3) F^ = argmin||[$ihi,...,$jhj] - G ||f s.t. supp(H) C 

4) = G - 

5) Stop if Ylj lb < e- 
Output: F^. 


Note that the linear constraint C is not enforced in Algorithm 2. The idea is to first find the 
support of the multi-vectors without taking into account of the linear constraint, and, in the second 
stage, follow the support recovery with least squares 

(38) F* = arg min ||G — [$ihi,..., $jhj]||F, s.t. supp(H) C supp(F°°), £H = 0 

H 
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where F°° is the output of Algorithm 2. 

For more discussion and applications of constrained joint sparsity, the reader is referred to 
Fannjiang 2013a where the performance guarantees similar to Theorem[^and Theoremj^are proved 
for constrained joint sparsity. 


5. Fresnel diffraction with point objects 


A major problem with discretizing the object domain shows up when the objects are point-like. 
In this case it is unrealistic to assume the objects are located exactly on the grid as the forceful 
matching between the point objects and the grid can create detrimental errors. Without additional 
prior information the gridding error due to the mismatch between the point object locations and 
the grid points can be as large as the data themselves, resulting in a low Signal-to-Noise Ratio 
(SNR). 


We shall call the grid spacing I given in (28) the Resolution Length (RL), which is the natural 
unit for resolution analysis. In the RL unit, the object domain grid becomes a subset of the integer 
grid Z. 

In the case of point objects, to refine the standard grid and reduce discretization error we consider 
a fractional grid 


(39) 


Z/F = {j/F : j G Z} 


where F G N is called the refinement factor. The random partial Fourier matrix (25) now takes the 
form 


(40) 


$ = 


^—i2TT^jk/F 


where fj G [—1/2,1/2] are independent uniform random variables. In the following numerical 


examples, we shall consider both deterministic (see (45)) as well as random sampling schemes. 

As shown in Fig. the relative gridding error ||d||/||$f|| is roughly inversely proportional to 
the refinement factor F. 

Fig. § shows the coherence pattern [//(j, k)] of a 100 x 4000 matrix (|40[) with F = 20 (left 


panel). The bright diagonal band represents a heightened correlation (pairwise coherence) between 
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relative gridding error versus the refinement tractor 



Figure 4. The relative gridding error is roughly inversely proportional to the re¬ 
finement factor. (Fannjiang and Liao 2012a. Copyright @2012 Society for Industrial 
and Applied Mathematics. Reprinted with permission. All rights reserved) 


pairwise coherence pattern 



100*4000 matrix with F = 20 & coherence = 0.99566 


coherence versus the radius of exoluded band 



Figure 5. Coherence pattern [ii{j,k)] for the 100 x 4000 matrix with F = 20 
(left). The off-diagonal elements tend to diminish as the row number increases. The 
coherence band near the diagonals, however, persists, and has the average profile 
shown on the right panel where the vertical axis is the pairwise coherence averaged 
over 100 independent trials and the horizontal axis is the distance between two 
object points (Fannjiang and Liao 2012a. Copyright @2012 Society for Industrial 
and Applied Mathematics. Reprinted with permission. All rights reserved). 


a column vector and its neighbors on both sides (about 30). The right panel of Figure shows 

a half cross section of the coherence band across two RL, averaged over 100 independent trials. 

In general sparse recovery with large F exceeds the capability of currently known algorithms as 

the condition number of the 100 x 30 submatrix corresponding to the coherence band in Figure 
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easily exceeds 10^®. The high condition number makes stable recovery impossible. While Figure 
is typical of the coherence pattern of one-dimensional sensing matrices, the coherence pattern 
for two or three dimensions is considerably more complicated depending on how the objects are 
vectorized. 


5.1. BLOOMP. To overcome the conundrum of highly coherent sensing matrix due to a refined 
grid, we have to go beyond the coherence parameter and study the coherence pattern of the sensing 
matrix. 

The coherence pattern of a sensing matrix can be described in terms of the notion of coherence 
band defined below. Let rj > 0. Define the ry-coherence band of the index k as 

(41) Brj{k) = {i \ fj,{i,k) > T]}, 
and the double coherence band as 

(42) Bf(A:) = B,{B,{k)) = Uj^B,(k)B,{j) 

The first technique for taking advantage of the prior information of well separated objects is 
called Band Exclusion (BE) and can be easily embedded in the greedy algorithm. Orthogonal 
Matching Pursuit (OMP). 

To imbed BE into OMP, we make the following change to the matching step 

Wx = argmini I, i ^ n = l,2,.... 

meaning that the double ?y-band of the estimated support in the previous iteration is avoided in the 
current search. This is natural if the sparsity pattern of the object is such that Br^{j),j G supp(f) 
are pairwise disjoint. We call the modified algorithm the Band-excluded Orthogonal Matching 
Pursuit (BOMP) as stated in Algorithm 3. 
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Algorithm 3. Band-Excluded Orthogonal Matching Pursuit (BOMP) 
Input: ^,g,T] > 0 

Initialization: = 0, = g and = 0 

Iteration: For j = 1,s 

1) ima.x = argmaxj | \,i ^ 

2) = 5^-1 U { i ^^} 

3) P = argminh ||$h — g ||2 s.t. supp(h) CS^ 

4) r-^ = g — 

Output: P. 


The following theorem gives a (pessimistic) performance guarantee for BOMP. 


Theorem 4. (Fannjiang and Liao 2012a) Let f be s-sparse. Let rj > 0 be fixed. Suppose that 


(43) 

Br,{i)FB)^){j)=%, 

\/i,j G supp{I) 

and that 



(44) 

ri{5s-A)^ + 

/min 

5||b||2 

2/min 

where 

/max = max|/fc|, 
k 

/min = min|/fc| 
k 


Let f® be the BOMP reconstruction. Then supp{i^) C Brj{supp{p) and moreover every nonzero 
component of f® is in the t]- coherence band of a unique nonzero component of f. 


Remark 3. Condition (fS) means that BOMP guarantees to resolve 3 RL. In practice, BOMP can 
resolve objects separated by close to 1 RL when the dynamic range is nearly 1. 


Remark 4. A main difference between Theorem and Theorem lies in the role played by the 


dynamic range /max//min and the separation condition (43). 

Another difference is approximate recovery of support in Theorem versus exact recovery of 

support in Theorem^ (a). In contrast to F-independent nature of approximate support recovery, 

exact support recovery would probably be highly sensitive to the refinement factor F. That is, as 

19 









F increases, the chance of missing some points in the support set also increases. As a result, the 
error of reconstruction ||/* — /II 2 tends to increase with F (as evident in Fig. [^. 

A main shortcoming with BOMP is in its failure to perform even when the dynamic range is even 
moderately greater than unity. To overcome this problem, we introduce the second technique: the 
Local Optimization (LO) which is a residual-reduction technique applied to the current estimate 
of the object support (Fannjiang and Liao 2012a). 


Algorithm 4. 

Local Optimization (LO) 

Input:$,g,7/ > 0,5° = {ii,... ,4}. 


Iteration: For j = 1, 2,..., k. 


1) P = arg minh ^•h — g 2 , supp(h) = 

(5^-^{4})u{z'},i' gb,({7,}). 

2) = supp(P). 


Output: 5^. 



In other words, given a support estimate , LO hne-tunes the support estimate by adjusting each 
element in within its coherence band in order to minimize the residual. The object amplitudes 
for the improved support estimate are obtained by solving the least squares problem. Because of 
the local nature of LO, the computation is efficient. 

Embedding LO in BOMP gives rise to the Band-excluded, Locally Optimized Orthogonal Match¬ 
ing Pursuit (BLOOMP). 


Algorithm 5. Band-excluded, Locally Optimized Orthogonal Matching Pursuit (BLOOMP) 
Input: g, 7 / > 0 

Initialization: = 0, = g and = 0 

Iteration: For j = 1,..., s 

1) Wx = argmaxj | \,i ^ 

2) = LO(5'^“^ U {imax}) where LO(5'^“^ U {imax}) is the output of Algorithm 4 
with U {i max } as input. 

3) O = argmiuh ||$h — g ||2 s.t. supp(h) G 

4) rJ = g - 
Output: P. 


20 








The same BLO technique can be used to enhance the other well known iterative schemes such 
as SP, CoSaMP (Needell and Tropp 2009), Compressed Iterative Hard Thresholding (IHT) (Blu- 
mensath and Davies 2009, Blumensath and Davies 2010)and the resulting algorithms are denoted 
by BLOSP, BLOCoSaMP and BLOIHT, respectively, in the numerical results below. We refer the 
reader to Fannjiang and Liao 2012a for the details and descriptions of these algorithms. 

MATLAB code of Algorithm 3.5 is available on-line at 

https://www.math.ucdavis.edu/~ fannjiang/home/codes/BLOOMPcode. 

5.2. Band-excluding thresholding. A related technique that can be used to enhance BPDN/Lasso 
for off-grid objects is called the the Band-excluding, Locally Optimized Thresholding (BLOT). 


Algorithm 6. Band-excluding, Locally Optimized Thresholding (BLOT) 

Input: f = (/i,..., /at), ^,g,V> 0- 
Initialization: = 0. 

Iteration: For j = 1, 2,..., s. 

1) ij = arg max|/fc|,A: 0 

2) = 5^-1 U {ij}. 

Output: P = argmin ||$h — g|| 2 , supp(h) C LO(S'®) where LO is the output of Algorithm 4. 


5.3. Numerical examples. For numerical demonstration in Fig. IE we use deterministic, equally 
spaced sampling with 


(45) 




and $ G £_MxFM _ 150, T = 50 to recover 20 randomly distributed and randomly phased 

point objects (spikes) separated by at least 4 RL. 

Fig. § (a)(b) show how the BLO technique corrects the error of OMP due to the unresolved 

grid. In particular, several misses are recaptured and false detections removed. Fig. (c) (d) show 

how the BLOT technique improves the BPDN estimate. In particular, BLOT has the effect of 

“trimming the bushes” and “growing the real trees”. Fig. [^a through c shows the relative error 

of reconstruction as a function of F by OMP, BPDN, BLOOMP and BPDN-BLOT with the same 
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OMP 


BLOOMP 



(a) OMP 



residual = 5.0332% L1 error = 158.9065% 0.1 -filtered error = 29.3316% 


(c) BPDN 



dis = 0.04 residual = 4.4366% LI error = 38.9837% 0.1 -filtered error = 6.0039% 

(b) BLOOMP 


BPDN-BLOT 



dis = 0.04 residual = 4.4366% LI error = 38.9837% 0.1-filtered error = 6.0039% 


(d) BP-BLOT 


Figure 6 . Reconstruction by (a) OMP, (b) BLOOMP, (c) BPDN and (d) BPDN- 
BLOT of the real part of 20 randomly phased spikes with F = 50, SNR = 20 
(Fannjiang and Liao 2012b. Reprinted with permission). 


set-up and three different SNRs. For all SNRs, BLOOMP and BPDN-BLOT produce drastically 
less errors compared to OMP and BPDN. 

The growth of relative error with F reflects the sensitivity of the reconstruction error alluded to 
in Remark Note that the reconstruction error in the discrete norm can not distinguish how far 
off the recovered support is from the true object support. The discrete norm treats any amount 
of support offset equally. An easy remedy to the injudicious treatment of support offset is to use 
instead the filtered error norm ||f^ — f^||, where f,j and are, respectively, f and f* convoluted with 
an approximate delta-function of width 2r]. 
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(a) SNR=100, ?7 = 0 


(b) SNR=20, 77 = 0 


(c) SNR=10, rj — 0 



(d) SNR=100, n = 0.05^ 


(e) SNR=20, 77 = 0.05£ 


(f) SNR=10, 77 = 0.05f 


Figure 7. Relative errors in reconstruction by OMP, BLOOMP, BP and BP-BLOT 
as F varies (top) without or (bottom) with filtering (Fannjiang and Liao 2012b. 
Reprinted with permission). 


Clearly the filtered error norm is more stable to support offset, especially if the offset is less 
than rj. If every spike of is within rj distance from a spike of f and if the amplitude differences 
are small, then the ry-filtered error is small. As shown in Fig. (d)(e)(f), averaging over rj = 5% 
RL produces acceptable filtered error for any refinement factor relative to the external noise. This 
suggests that both BPDN-BLOT and BLOOMP recover the object support on average within 5% 
of 1 RL, a significant improvement over the theoretical guarantee of Theorem]^ 


Next we consider the unresolved partial Fourier matrix (40) with random sampling points to 

demonstrate the flexibility of the techniques. Let G [—1/2,1/2], j = 1,...,M be independent 

uniform random variables with M = 100, N = 4000 and F = 20. The test objects are 10 randomly 

phased and distributed objects, separated by at least 3 RL. As in Theorem]^ a recovery is counted 

as a success if every reconstructed object is within 1 RL of the object support. 
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success probability versus SNR when dynamic range = 1 



success probability verus dynamic range when SNR = 33 



30 40 50 

dynamic range 


Figure 8. Success probability versus (left) SNR for dynamic range 1 and (right) 
dynamic range for SNR = 33. Here LOOMP is a simplified version of BLOOMP 
and has nearly identical performance curves (Fannjiang and Liao 2012a. Copyright 
@2012 Society for Industrial and Applied Mathematics. Reprinted with permission. 
All rights reserved). 


Fig. compares the success rates (averaged over 200 trials) of the BLO-enhanced schemes 
(BLOOMP, BLOSP, BLOCoSaMP, BLOIHT) and BLOT-enhanced scheme (Lasso-BLOT). Lasso- 


BLOT is implemented with the regularization parameter 


(46) 


A = 0.5@log N (black curves with diamonds) 


or 


(47) 


A = @ 2 log N 


(black curves with stars) 


(Chen et al. 2001). The empirically optimal choice (46) (labelled 


much improved performance over the choice (47). Clearly, BLOOMP 
stability and dynamic range among all tested algorithms. 


as Lasso-BLOT (0.5)) has a 
is the best performer in noise 


5.4. Highly redundant dictionaries. Our discussion in Section so far is limited to point-like 

objects. But the methods presented above are also applicable to a wide variety of cases where the 

objects have sparse representations by redundant dictionaries, instead of orthogonal bases. 

Suppose that the object is sparse in a highly redundant dictionary, which by definition, tends to 

represent an object by fewer number of elements than a non-redundant one does. For example, one 

can combine different orthogonal bases into a dictionary that can sparsify a wider class of objects 
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Figure 9. Single-pixel camera block diagram (http://www.dsp.ece.rice.edu/cscainera/) 


than any individual base can. On the other hand, a redundant dictionary tends to produce a larger 
coherence parameter and be ill suited for compressive sensing. This is the same kind of conundrum 
about off-grid point-like objects. 

One of the most celebrated examples of optical compressive sensing is the Single-Pixel Camera 
(SPC) depicted in Fig. In SPC, measurement diversity comes entirely from the Digital Mi¬ 
cromirror Device (DMD) instead of sensor array. The DMD consists of an array of electrostatically 
actuated micro-mirrors. Each mirror can be positioned in one of two states (±12°). Light reflected 
from mirrors in the ±12°-state only is then collected and focused by the lens and subsequently 
detected by a single optical sensor. For each and every measurement, the DMD is randomly and 
independently reconfigured. The resulting measurement matrix A has independently and identi¬ 
cally distributed entries. 

Suppose that the object is sparse in terms of a highly redundant dictionary. For simplicity of 
presentation, consider an ID object sparse in an over-complete Fourier frame (i.e. a dictionary that 
satishes the frame bounds Daubechies (1992) ) with entries 

(48) fe = l,. j = l,...,RF, 

that includes harmonic as well as non-harmonic modes as its columns, where F is the redundant 
factor and R is a large integer. In other words, the object can be written as 'J'f with a sufficiently 
sparse vector f. The final sensing matrix then becomes 


$ = A^. 


(49) 


25 
















coherence versus the radius of excluded band coherence versus the radius of excluded band 




Figure 10. The coherence bands of the redundant Fourier frame 'I' (left) and 
$ = (right), the latter being averaged over 100 realizations of A (Fannjiang 
and Liao 2012a. Copyright (6)2012 Society for Industrial and Applied Mathematics. 
Reprinted with permission. All rights reserved). 


The coherence bands of 'I' and $ are shown in Figure [l0| from which we see that like Fig. [^the 
coherence radius is less than 1 RL. The same BLO- and BLOT-based techniques can be applied to 


(49), see Fannjiang and Liao 2012a for numerical results and performance comparison with other 
techniques for off-grid objects (Candes et al. 2011, Candes and Fernandez-Granda 2013, Candes 
and Fernandez-Granda 2014, Duarte and Baraniuk 2013, Tang et al. 2013). 


6. Fresnel diffraction with Littlewood-Paley basis 

Opposite to the localized pixel basis, the Littlewood-Paley basis is slowly decaying, nonlocal 
modes based on the wavelet function 


(50) 


'6(x) = (nx) ^(sin (27ra:) — sin (vrx)) 


which has a compactly supported Fourier transform 


(51) 

The following functions 


I 0, otherwise. 


'il>p,g{x) = 2 - q), 
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(52) 


p,q eZ 













form an orthonormal wavelet basis in L^(M) (Daubechies 1992). Expanding the masked object V 
(20) in the Littlewood-Paley basis we write 


(53) V(x)= ^ Vp^glpp^qix). 

p,q&Z 

The main point of the subsequent discussion is to design a sampling scheme such that the 
resulting sensing matrix has desirable compressive sensing properties (Fannjiang 2009). 

Let {2P : p = —p*, —p* + 1, ...,p*} be the dyadic scales present in ( [5^ , {q : |g| < A^p} the modes 
present on the scale 2^ and 2Mp + 1 the number of measurements corresponding to the scale 2^. 
Let 


p'-i 


(54) 


k = + 1) + W\ ^ \p\ 


< p* 


J=-p* 


be the index for the sampling points. Throughout this section, k is determined by p',q' by (54). 
Let Xfc be the sampling points and set the normalized coordinates 


(55) 


XkU}£ 

2ttzo 


= ik, k = 1, ...M 


where, as shown below, £ is a resolution length and £ [~l/2jl/2] are determined below, c.f. 


(23). This means that the aperture (i.e. the sampling range of Xk) is again given by (27). 
Let g = (pk) be the data vector with 


9k 




Direct calculation with ([^ and (55) then gives 


(56) gk = ^ 2?’/Vp,,e-'2<'=^“^2^V(ar'2?’), fc = l,...,M. 

p,q& 

Let f = {fi) be the object vector with 


fi = {-iy2P/%,q 


where the indices are related by 


p-i 

I = ^ {2Nj + 1) + g. 

3=-p* 
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Suppose that 


(57) 


£<2-p.-i 


i.e. 2i is less than or equal to the smallest scale in the wavelet presentation (53) 
Let Cp',q' be independent, uniform random variables on [—1/2,1/2] and let 


(58) 


^ £ J l/2 + Cp',g'5 Cp',q' ^ [0)1/2] 

2?>' I + e [_i/2,0] 


where k is determined by (54). By the assumption (57), we have 


G [-l/2,l/2], Vp'>-p,. 


More specifically, by (55), we have 


Xk G 


27rzo 

u;2P' 



1' 


T 

( 

~2 

U 

2’ ^ 


i.e. the sampling regions for different dyadic scales indexed by p' are disjoint with the ones for 
the smaller scales on the outer skirt of the aperture, taking up a bigger portion of the aperture. 
The resulting sampling points are geometrically concentrated near (but not exactly at) the center 
of the aperture. 

Let the sensing matrix elements be 


(59) 


= (-l)''^(a2^ri)e-'2-«^W. 


(60) 


We claim that = 0 for p ^ p'. This is evident from (58) and the following calculation 

= 2P-P' 


1/2 + Cp',q'i Cp',q' G [0, 1/2] 
“1/2 + Cp',ij') Cp',q' G [—1/2,0]. 


For p ^ p' the absolute value of (60) is either greater than 1 or less than 1/2 and hence (60) is 
outside the support of il) . 


On the other hand, ioi p = p', (60) is inside the support of ifj and so 
(61) ^k,l = Jg'l < Mp, \q\ < Np 
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which constitute the same random partial Fourier matrix that we have seen above. In other words, 


under the assumption (57) the sensing matrix $ = [$*;,/] e with N = Y1 p|<p, (2iVp + 1) and 

M = X]|p|<p^(2Mp + 1), is block-diagonal with each block (indexed by p) in the form of random 
partial Fourier matrix, representing the sensing matrix on the dyadic scale 2^. 


7. Near-field diffraction with Fourier basis 


Consider near-held diffraction by a periodic, extended object (e.g. diffraction grating) where the 
evanescent modes as well as the propagation modes are taken into account. Since we can not apply 
the paraxial approximation, we resort to the Lippmann-Schwinger equation (§. 

Suppose the masked object function is sparse in the the Fourier basis 

OO 

(62) V{x)= Y, 

j=-oo 

where L is the period and only s modes have nonzero amplitudes. Suppose that Vj = Q for 

The 2D Green function can be expressed by the Sommerfeld integral formula 


(63) 
where 

(64) 




^ii^(\z\l3{a)+xa) 


da 


r = {z, x) 


/3(a) = 




a^ 


|a| < 1 


i\/a^ — 1, |a| > 1 


(Born and Wolf 1999). The integrand in (63) with real-valued /3 (i.e. |a| < 1) corresponds to the 
homogeneous wave and that with imaginary-valued /3 (i.e. |a| > 1) corresponds to the evanescent 
(inhomogeneous) wave which has an exponential-decay factor Likewise the 3D Green 

function can be represented by the Weyl integral formula (Born and Wolf 1999). 

The signal arriving at the sensor located at (0, x) is given by the Lippmann-Schwinger equation 


with (63) 
(65) 




G{zo,x-x')V{x')dx' = 
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where 


( 66 ) 


CTi = 


27rj 
Luj ’ 


/3j = /3iaj). 


The subwavelength structure is encoded in Vj with aj > 1 corresponding to the evanescent modes. 

Let {0,Xk),Xk = ^kL,k = be the coordinates of the sampling points where £ 

[—1/2,1/2]. In other words, L is also the aperture (i.e. the sampling range for x^)- To set 
the problem in the framework of compressed sensing we set the vector f = (fj) G as 


(67) 


■ iuzoPj 

T_ V- 

2uj/3j 


To avoid a vanishing denominator in ( [67^ , we assume that aj ^ 1 and hence 7 ^ 0,Vj G Z. This 
is the case, for instance, when Loj/{2'k) is irrational. 

This gives rise to the sensing matrix $ with the entries 


( 68 ) 


^kj 






which again is the random partial Fourier matrix. 


A source of instability lurks in the expression (67) where f3j may be complex-valued, corre¬ 


sponding to the evanescent modes. Stability in inverting the relationship (67) requires limiting the 


number of the evanescent modes involved in (67). Here the transition is not clear-cut, however. 
For example, if we demand that 


(69) 


JuizoP^ 


I > e 


-27r 


as the criterion for stable modes, then the stable modes include \aj\ < 1 as well as \aj\ > 1 such 
that 

(70) u}\fij\zQ < 27r 


or equivalently 
(71) 





1 
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In other words, the number of stably resolvable modes is proportional to the probe frequency and 
inversely proportional to the the distance zq between the sensor array and the object. As zq drops 
below the wavelength, the subwavelength Fourier modes of the object can be stably recovered. This 
is the idea behind the near-held imaging systems such as the scanning microscopy. 


8. Inverse scattering 

In the inverse scattering theory, the scattering amplitude is the observable data and the main 
objective then is to reconstruct ly from the knowledge of the scattering amplitude. 

8.1. Pixel basis. To obtain a sensing matrix with compressive sensing properties, we hrst make 
the Born approximation in Q and neglect the scattered held u® on the right hand side of Q. 
Our purpose here is to demonstrate how to coordinate the incidence direction and the sampling 
direction and create a favorable sensing matrix. 

Consider the incidence held 


(72) 


tt^r) = e 




where d is the incident direction. Under the Born approximation, we have from Q that 
(73) 


u 


A(f, d) = A{s) = ^ / i/(r')e-*‘^’' '"dr' 

dvr ./rad 


where s = r — d is the scattering vector. 


We proceed to discretize the continuous system (73) as before. Consider the discrete approxi¬ 
mation of the extended object u 


(74) 


where 


(75) 


Mr) = Kj- q)i^(^q) 


5(r) = 


1, r e [-2, 


1 112 
2’ 2J 


0, else. 


is the pixel basis. 

Dehne the target vector f = (fj) G with fj = z^(.^p),p = (pi,P 2 ) £ = (Pi ~ 1)V^ + P 2 - 

Let oji and d/ be the probe frequencies and directions, respectively, and let r; be the sampling 
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directions for I = 1,M. Let g be the data vector with 

47r^(r; - df) 

Then the sensing matrix takes the form 

(76) ci={qi,q2)€Z%, j = (q, - 1)^ + q^. 

8.2. Sampling schemes. Our strategy is to construct a sensing matrix analogous to the random 
partial Fourier matrix. To this end, we write the (Z,j)-entry of the sensing matrix in the form 


where ^i,(i are independently and uniformly distributed in [—1,1]. Write in the polar 

coordinates pi , 4>i as 

(77) {(iXi) = Pi{cos4>i,sm(pi), pi = \J if + Cz^ < V2 

and set 


uji{cos6i — cos 6i) = V2piQcos(f)i 

uji{sm9i — sin 6i) = V^piQsincpi 


where is a parameter to be determined later (91). Equivalently we have 


(78) 

— V2uji sin 

(79) 

y/2'^l sin 


cos 


= Qpi sin 4>i. 


2 2 

This set of equations determines the single-input-(0;, cu;)-single-output-0; mode of sampling. 


The following implementation of (78)-(79) is natural. Let the sampling angle 9i be related to the 
incident angle 9i via 


(80) 


9i + 9i = 2(t)i + TT, 
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and set the frequency oji to be 


(81) 


ui = 


Qpi 


V2 


sin 


ei-ei 


Then the entries (76) of the sensing matrix $ have the form 


(82) 




By the square-symmetry of the problem, it is clear that the relation (80) can be generalized to 


(83) 


9i + 9i = 2(t)i -7 r/TT, 7/ E Z. 


On the other hand, the symmetry of the square lattice should not play a significant role and hence 


we expect the result to be insensitive to any fixed ?? E M, independent of I, as long as (81) holds. 
Indeed this is confirmed by numerical simulations. 

Let us focus on two specific measurement schemes. 


Backward sampling. This scheme employs 11—band limited probes, i.e. w; E [—11,0]. This and (81) 
lead to the constraint: 


(84) 


sm • 


01-9i 


> 


Pi 

V2- 


The simplest way to satisfy (80) and (84) is to set 


(85) 

( 86 ) 


fil = 91 = 01 + TT, 

Qpi 


UJl = 


^/2 


I = l,...,n. In this case the scattering amplitude is always sampled in the back-scattering direc¬ 
tion. This resembles the synthetic aperture imaging which has been previously analyzed under the 
paraxial approximation in Fannjiang et al. 2010. In contrast, the forward scattering direction with 


9i = 9i almost surely violates the constraint (84). 


Forward sampling. This scheme employs single frequency probes no less than 11: 


(87) 


ui = 7 II, 7 > 1 , 1 = 1, ...,n. 
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To satisfy (83) and (81) we set 


TI7T 

( 88 ) 0i = (pi + — + arcsin 

(89) ei = (pi + ^ - arcsin 

with r] € Tj. The difference between the incident angle and the sampling angle is 

(90) 6i — 9i = 2 arcsin 

which diminishes as 7 —)> oo. In other words, in the high frequency limit, the sampling angle 
approaches the incident angle. This resembles the setting of the X-ray tomography. 

In summary, let be independently and uniformly distributed in [—1,1] and let {pi,<pi) be 
the polar coordinates of (61 C/); be. 


Pi 

7V2 


Pi 

-fV2 


{Cl, Cl) = Pi{cos(pi,sm(pi). 


Then with with 

(91) M = TrlV2 

both forward and backward samplings give rise to the random partial Fourier sensing matrix. 


8.3. Coherence bounds for single frequency. As in Section we let the point scatterers 
be continuously distributed over a finite domain, not necessarily on a grid. Any computational 
imaging would involve some underlying, however refined, grid. Hence let us assume that there is 
an underlying, possibly highly refined and unresolved, grid of spacing i <C (the reciprocal of 
probe frequency). 

We shall focus on the monochromatic case with cvi = uj,l = 1,M. 


Recall the sensing matrix continues of the form (76) which now becomes 


(92) 


iPij = j = (p^_i)ViV + p2, PGZ 


N- 


In other words, the measurement diversity comes entirely from the variations of the incidence and 

detection directions. We assume that the n incident directions and the m detection directions are 
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each independently chosen according to some distributions with the total number of data M = nm 
hxed. 


Theorem 5. (2D case). Suppose the incident and sampling angles are randomly, independently 
and identically distributed according to the probability density functions ffO) G and P{9) G , 
respectively. Suppose 


(93) 


e,K>0. 


Set L = — q| for any p,q G Then the sensing matrix satisfies the pairwise coherence 

bound 


(94) 


( i V2k\ ( 3 y/2K 


with probability greater than (1 — e)^ where 


(95) 

(96) 

with a positive constant c. 


p} < c{l + uL) supg {\f{e)\,\^f{e)\} , 
< c(l + u;L)-V 2 ^ 


In 3D, the coherence bound can be improved with a faster decay rate in terms of cvL S> 1 as 
stated below. 


Theorem 6. (3D case). Assume (93). Suppose the incidence and sampling directions, parametrized 
by the polar angle 6 G [0,7r] and the azimuthal angle (p G [0,27r], are randomly, independently and 
identically distributed. Let f\9) G and f^{0) G be the marginal density functions of the 
incident and sampling polar angles, respectively. 

Let L = £|p — q|. Then the sensing matrix satisfies the pairwise coherence bound 


(97) 


^p.q< (p}+^^^ (p- + 

\ ^/n I \ y/m I 


V2k \ 


with probability greater than (1 — e)^ where 


(98) 

(99) 


p^ < c{l + u}L) ^supe {|/‘(6')|, |^/'(6')|} 

p^<c{l + uL)-^snvg{\nB)\,\l,r{e)\}. 
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Figure 11. Two instances of BOMP reconstruction: red circles are the exact loca¬ 
tions, blue asterisks are recovered locations and the yellow patches are the coherence 
bands around the objects. 

Remark 5. The original statements of the theorems (Fannjiang 201 Oh, Theorems 1 and 6) have 
been adapted to the present eontext of off-grid objects. The original proofs, however, carry over 
here verbatim upon minor change of notation. 

Remark 6. When the sampling directions are randomized and the incidence directions are deter¬ 


ministic, then the coherence bounds (94) and (97) hold with the first factor on the right hand side 
removed. 


( 100 ) 

( 101 ) 


According to Remark]^ we have the pairwise coherence bound: 

d 


(2D) /ip,q <c(l + a;L) sup 11/®(6') 

e I 

(3D) /ip,q <c{l + LvL)~^supi.\f{e)\, 

e { 


dO 


no) 


+ 


0 


no) 


+ 


V2K 

y/M 

V2K 

y/M 


which is an estimate of the coherence pattern of the sensing matrix. Hence, if L is unresolvable (i.e. 
(wL < 1), the corresponding pairwise coherence parameter is high and when if L is well-resolved 
(i.e. cjL ^ 1) the corresponding pairwise coherence parameter is low. A typical coherence band 


has a coherence radius 0{u) ^) according to (lOO)-(lOl). 


Therefore, if the point objects are well separated in the sense that any pair of objects are larger 

than u)~^ then the same BLO- and BLOT-based techniques discussed in Section]^ can be used to 

recover the masked object support and amplitudes. For a simple illustration, Figure pT] shows two 

instances of reconstruction by BOMP. The recovered objects (blue asterisks) are close to the true 

objects (red circles) well within the coherence bands (yellow patches). 
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9. Inverse multiple scattering 


In this section, we present an approach to compressive imaging of multiply scattering point 
scatterers. First consider the multiple scattering effect with just a single illumination, i.e. ra = I 
and M = m. 

Note that the original object support is the same as the masked object support. With the support 
accurately recovered, let us consider how to unmask the objects and recover the true objects. 

Define the incidence and full field vectors at the locations of the objects: 

u‘ = (tt‘(ri), ...,u‘(rs))^ G C® 

u = (u(ri), ...,R(r^))'^ G 


Let r be the s x s matrix 


and V the diagonal matrix 




V = diag(i/i,...,Us). 

The full field is determined by the Foldy-Lax equation (Mishchenko et al. 2006) 


( 102 ) 


u = u' + w^FVu 


from which we obtain the full field 


(103) 


and the masked objects 


(104) 


u 


= (I - oi^rv) u‘ 


f = vu = v(i-w2rv) 

= (l-a;2vr)“Vu' 


provided that w ^ is not an eigenvalue of FV. 


Hence by (104) we have 


(105) 


(I - w^vr) f = vu\ 
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The true objects v can then be recovered by solving (105) as 


(106) 


f 

a;2rf + ui 


where the division is carried out entry-wise (Hadamard product). 


9.1. Joint sparsity. With the total number of data M = nrti fixed the coherence bounds (94) and 


(97) is optimized with n ^ m ^ \fM. To take advantage of this result, we should deploy multiple 


incidence fields for which the formula (106) is no longer valid 


Multiple illuminations give rise to multiple data vectors gj and multiple masked object vectors 
fj, j = 1, ...,n each of which is masked by a unknown field Uj. However, all masked object vectors 
give rise to the same sensing matrix 


J = (pi-l)v/iv+p2, PGZ; 


'N- 


Since every masked object vector shares the same support as the true object vector, this is a 


suitable setting for the application of joint sparsity techniques discussed in Sections 4^ and 4.3 
Compiling the masked object vectors as F = [fi,. ..,£„] G ([^mxn data vectors as G = 

[gi) ••■jgn] £ we obtain the imaging equations 


(107) 


G = + E 


where E accounts for noise. When the true objects are widely separated, we have two ways to 
proceed as follows. 


1) BPDN-BLOT for joint sparsity. In the first approach, we nse BPDN for joint sparsity (37) 
with = 0 to solve the imaging equation (107). Let F* = (fi*, ...,fn*) be the solution. 


We then apply the BLOT technique (Algorithm 5) to improve F*. In order to enforce the joint 

sparsity structure, we modify Algorithm 5 as follows. 

First, we modify the LO algorithm to account for joint sparsity. 
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Algorithm 7. LO for joint sparsity 
Input: > 0, 5° = 

Iteration: For k = 1,2,..., s. 

1) = argmin||[$ihi,...,$nh„] - G||f s.t. UjSupp(hj)C(5^“^\{zfc}) U {^j, 4 £ Br,{{ik})- 

2) = supp(F*^). 

Output: S^. 


Next, we modify the BLOT algorithm to account for joint sparsity. 


Algorithm 8 . BLOT for joint sparsity 
Input: ^ 1 , ..,^n,G,rj > 0. 

Initialization: = ^. 

Iteration: For k = 1,2,..., s. 

1) ik = arg maxj ||fjH 2 , k 0 

2) 5^ = 5^-1 U {4}. 

Output: F* = argmin ||[$ihi,$nhn] — G||f, UjSupp(hj)CJLO(5^) where JLO(5^) 
is the output of Algorithm 7 with the s-th iterate of BLOT as input. 


2) BLOOMP for joint sparsity. In the second approach, we propose the following joint sparsity 
version of BLOOMP. 


Algorithm 9. BLOOMP for joint sparsity 

Input: $i,...,$„,G,?7 > 0 
Initialization: F® = 0, = G and 5 ® = 0 

Iteration: For k = 1,..., s 

1) imax = argmaxj^^^^ |<I>]^ B^\s^~^), where <I)j • = conjugate transpose of coli(^i)- 

2) = JLO(5^“^ U {imax}) where JLO is the output of Algorithm 7. 

3) [f4...,fn] = argmiuH ||[^ihi,...,$„h„] - G||f s.t. Ujsupp(hj) 

4) [r4...,r^] = G-[$iff,...,$„f„^] 

Output: F* = [ff,...,f^]. 
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After the first stage of either approach, we obtain an estimate of the object support as well as 
the amplitudes of masked objects. In the second stage, we estimate the true object amplitudes. If 


we use the formula (106) for each incident wave ul. 


we end up with n amplitude estimates 


‘■j* 


+ ui 


j = 1, ...,n 


that are typically inconsistent. Least squares is the natural way to solve this over-determined 
system and obtain the object estimate 


= argmin^ ||(a;^rfj* -|- u‘ )v — 
i=i 

10. Inverse Scattering with Zernike basis 


In this section, we discuss a basis for representing extended objects in the scattering geometry 
and its application to compressive inverse scattering. We shall make the Born approximation. 

A well known orthogonal basis for representing an extended object with a compactly support 
(e.g. the unit disk) is the product of Zernike polynomials i?™ and trigonometric functions 


(108) = K"(pcos0,psin0) = x ^ + y‘^<l 


where m G Z,n G N, n > \m\ and n — |m| is even. We refer to as the Zernike functions of 
order (m,n) (Born and Wolf 1999). These Zernike functions are very useful in optics because the 
lowest few terms of a Zernike expansion have a simple optical interpretation (Dai and Mahajan 
2008). In addition, a Zernike expansion usually has a superior rate of convergence (hence sparser) 
compared with other expansions such as a Bessel-Fourier or Chebyshev-Fourier expansion (Boyd 
and Yu 2011 and Boyd and Petschek 2014). 

We show now that the Zernike basis also results in a better coherence parameter (hence better 
resolution) than the pixel basis. The Zernike polynomials are given explicitly by the formula 


(109) 


Kip) = 


[d{p 


r d 1 

n —|Tn| 

2 

[d{p^)\ 

- 


o. n+\m\ , „ , , n-|m| 

{p) 2 (P - 1) 2 
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which are n-th degree polynimials in p and normalized such that = 1 for all permissible 

values of m, n. The Zernike polynomials satisfy the following properties 




( 110 ) 

( 111 ) 




/ K{p)R::,{p)pdp = 

Jo + J-j 

f K{p)JUup)pdp = 

/o u 


where Jn+i is the (n + l)-order Bessel function of the first kind. As a consequence of (110), the 
Zernike functions satisfy the orthogonality property 


( 112 ) 


/ V^{x, y)V^'{x, y)dxdy = —^(5w(^nn'- 


Writing s = s{cos4>,sm<p), let us compute the matrix element for the scattering amplitude (73) 
as follows. 


(113) 


V;^{x,y)e-^‘^^<^’y^dxdy = 


1 /■27r 


' x^+y^<l 



10 JO 

U r2w 


^iujspcos (4>+e)ji^^py-imedgpdp 


0 Jo 


= 27ri"'e*”*‘^ / Jm{jJsp)Rn{p)pdp 

Jo 


by the definition of Bessel function 




1 r 

TTi^ Jq 


Az cos 6 


cos {m6)d6. 


Using the property (111), we then obtain from (113) that 


(114) [ = 27rz'"(-l)'^e*™^ 

J x^+y'^<l 


Jmd>'^n+lAs) 


(XS 


which are the sensing matrix elements with all permissible m, n. Note that the columns of the 
sensing matrix are indexed by the permissible m G Z, n G N with the constraint that n > |m| and 


n — m IS even. 


Let the scattering vector s = r — d be parametrized as 


Sjk = Sj(cos(/)A:,sin(/)fc), j,k = 1, ...,Vm 
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such that {4>k} are independently and identically distributed uniform random variables on [0,27r] 
and {sj} are independently distributed on [0, 2] according to the linear density function /(r) = r/2. 
As a result, Zj = ujSj are independently and identically distributed on [0, 2a;] according to a linear 
density function. 

Calculation of the coherence parameter between the columns corresponding to (m, n) ^ (m', n') 
gives the following expression 

J_ ^ Jn+l{^Sj) Jn/+i(a;Sj) \ ( 1 ^ 

ijjSj I I - 




UJSn 




Recall that for p,q G N 


(115) 


" U’ 


(Abramowitz and Stegun 1972, formula 11.4.6). For M 3> 1, we have by the law of large numbers 


(116) 


and 


(117) 


1 (a;Sj) 


E 


y/M ^ WSj 


E 


0084 


Jn+ljl-zr) Jn'+l{ujr) 
ur ojr 

r2w 


1 7^*^ dz 


1 ^27r 


= (5r, 


When m ^ m', the two columns are orthogonal and the pairwise coherence parameter is zero. 


When n ^ n', the right hand side of (116) becomes 0{u} ^) in view of (115) and the fact that 


the Bessel functions Jn{z) decay like z for z ^ 1. From (115) and (116) with n = n', we 


see 


that the 2-norm of the columns is 0(a;“^). After dividing (116) with n ^ n' hy the 2-norm of the 
columns the coherence parameter scales at worst like uj~^ (for m = m',n ^ n'). 

Notice that this decay date of the coherence parameter is faster than the behavior in 


(94)-(95). Hence, imaging with the Zernike basis possess better resolution capability than with the 


pixel basis, all else being equal. 
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11. Interferometry with incoherent sources 


In this last section, we discuss the compressive sensing application to optical interferometry in 


astronomy which has a similar mathematical structure to that of the inverse scattering (92) under 
the Born approximation. 

In astronomy, interferometry often deals with signals emitted from incoherent sources. In this 
section, we present compressive sensing approach to such a problem. With the help of the van 
Cittert-Zernike theorem, the sensing matrix has a structure not unlike what we discuss above. 

Suppose the field of view is small enough to be identified with a planar patch of the celestial 
sphere V, called the object plane. Let /(s) be the radiation intensity from the point s on the object 
plane V. Let n antennas be located in a square of size L on the sensor plane parallel to V with 
locations Lrj,j = 1, ...,n where rj G [0,1]^. Then by van Cittert-Zernike theorem (Born and Wolf 
1999) the measured visibility v{rj — r^) is given by the Fourier integral 


(118) 


{rj-rk)= f 

Jr 


Consider the discrete approximation of the extended object I with the pixel basis on the grid 


iZ% 


(119) 


h{r) = kI- q)-^(^q) 


where b is given in (75) and 


( 120 ) 


= {p = (Pl,P2) : P1,P2 = 1, ViV}. 


Substituting (119) into (118) we obtain the discrete sum 

N 


( 121 ) 


/ujiL 


v{rj - Tk) = fb f - rj) ) ^ 


'■{rj-rk)iL 


1=1 


where I, p are related by I = {pi — + P 2 and 

^ sin (vrO sin (vrr?) 

Vr^ VTT/ 

For every pair (j, k) of sensors we measure and collect the interferometric datum v{rj — r^) and we 
want to determine I from the collection of n(n — 1) real-valued data. 
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Let us rewrite eq. (121) in the form (10). In contrast to (28), we set 


( 122 ) 


TT 

ujL 


to account for the “two-way” structure in the imaging equation (121). Note that ^ is the resolution 
length on the celestial sphere and hence dimensionless. 

Let f = (/i) G be the unknown object vector, i.e. fi = Let g = {gi) G R^,M = 

n(n — l)/2, 


91 = 


SR [u(rj - rfc)] , I = {2n-j){j-l)/2 + k, j<fc = l,...,n 

& ((ffc — rj,)/2) I 3= [7;(rj — rfc)] , / = n(n — l)/2-|-(2n — j)(j — l)/2-|-A:, j<k = l,...,n 


be the data vector where 3^ and O' stand for, respectively, the real and imaginary parts. The sensing 
matrix $ G now takes the form 


(123) <l>ii = 


cos [27rp; • (r^ - r^)] 


i = (2n- j){j -l)/2 + k, j<k 


sin [27rp/• (rj - Tfc)], i = n{n - l)/2 + {2n - j){j - l)/2 + k, j<k 


which is no longer the simple random partial Fourier matrix for 2D as the baselines r,- — are 


related to one another. Nevertheless (123) has a similar structure to that of the inverse scattering 
(92) when the transmitters and receivers are co-located. Note that as (r^ — Yj)/2 G [—1/2,1/2]^ 


the denominator b ((r^ — rj)/2) in the definition of gi does not vanish. 

Next we give an upper bound for the coherence parameter. For the pairwise coherence for 
columns i, i' corresponding to p, p G we have the following calculation 

2 




n{n — 1) 


cos [27rp • {Tj - Tfc)] cos [27rp' • (r^ - r*,)] 


j<k 


First we claim: 


n(n — 1) 

1 

n(n — 1) 


/x(f,f') = 


-hsin [27rp • (r^ - r^)] sin [27rp' • (r^ - r^)] 


cos [27r(p - p') • {rj - rfc)] 

j<k 

^ COS [27r(p - p') • (rj - rfc)] 


_^_lly 

nin — 1) 11 

.7 = 1 




— n 
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This follows from the calculation 


g*27r(p-p')-rj 

j=l 


g*27r(p-p')-(rj-ri.) 

j¥=k 

^ cos [27r(p - p') • (rj - r^)] + i sin [27r(p - p) • (r^ - r^)] 

^ cos [27r(p - p') • {vj - Ffc)] 
j¥=k 


Some modification of the arguments for Theorems and [pleads to the following coherence bound. 


Theorem 7. Assume that the total number of grid point N satisfies the bound 

(124) N < 

with some constants 6 and K. Suppose that the sensor locations rj,j = are independent 

uniform random variables on [0,1]^. Then the coherence parameter pL satisfies the bound 

\2K‘^ — M 

(125) /r($) < '-^ 

n — 1 

with probability greater than 1 — 2e. 


In other words, with high probability the coherence parameter for the uniform distribution decays 
as n~^. A central problem in interferometry is the design of an optimal array, see Fannjiang 2013b 
for a discussion from the perspective of compressed sensing. 

Acknowledgements. Research is supported in part by US NSF grant DMS-1413373 and Simons 
Foundation grant 275037. 


References 

[1] Abramowitz, M. and I. Stegun. Handbook of Mathematical Functions (New York: Dover, New York, 1972). 

[2] Benedek, P. and R. Panzone, “The space H with mixed norm,” Duke Math. J. 28 (1961): 301-324. 

[3] Blumensath, T. and M.E. Davies, “Iterative hard thresholding for compressed sensing,” Appl. Comput. Harmon. 
Anal. 27 (2009): 265-274. 

[4] Blumensath, T. and M.E. Davies, “Normalized iterative hard thresholding: guaranteed stability and perfor¬ 
mance,” IEEE J. Sel. Top. Sign. Proc. 4 (2010): 298-309. 

[5] Born, M. and E. Wolf. Principles of Optics, 7-th edition. (Cambridge: Cambridge University Press, 1999). 

45 




[6] Boyd, J.P. and F. Yu, “Comparing six spectral methods for interpolation and the Poisson equation in a disk: 
radial basis functions, Logan-Shepp ridge polynomials, Fourier-Bessel, Fourier-Chebyshev, Zernike polynomials, 
and double Chebyshev series,” J. Comput. Phys. 230 (2011): 1408-1438. 

[7] Boyd, J.P. and R. Petschek, “The relationships between Chebyshev, Legendre and Jacobi polynomials: The 
generic superiority of Chebyshev polynomials and three important exceptions,” J. Sci. Comput. 59 (2014):l-27. 

[8] Boyd, S. and L. Vandenberghe. Convex Optimization. (Cambridge: Cambridge University Press, 2004). 

[9] Bruckstein, A.M., D.L. Donoho and M. Elad, “From sparse solutions of systems of equations to sparse modeling 
of signals,” SIAM Rev. 51 (2009): 34-81. 

[10] Candes, E. J., “The restricted isometry property and its implications for compressed sensing,” Compte Rendus 
de I’Aeademie des Seiences, Paris, Serie I. 346 (2008): 589-592. 

[11] Candes, E. J., Y.C. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant 
dictionaries,” Appl. Comput. Harmon. Anal.^l (2011): 59-73. 

[12] Candes, E. J., and C. Fernandez-Granda,“Super-resolution from noisy data” , Journal of Fourier Analysis and 
Applications 19 ( 6 ) (2013): 1229-1254. 

[13] Candes, E. J., and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” Commun. Pure 
and Applied Mathematics 67 ( 6 ) (2014): 906-956. 

[14] Candes, E. J. and T. Tao, “ Decoding by linear programming,” IEEE Trans. Inform. Theory 51 (2005): 4203- 
4215. 

[15] Candes, E. J., J. Romberg and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly 
incomplete frequency information,” IEEE Trans. Inform. Theory 52 (2006): 489-509. 

[16] Chambolle, A., “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20 
(2004): 89-97. 

[17] Chambolle, A. and P.-L. Lions, ’’Image recovery via total variation minimization and related problems, ” Numer. 
Math. 76 (1997): 167-188. 

[18] Chen, J. and X. Huo, “Theoretical results on sparse representations of mulitpie-measurement vectors,” IEEE 
Trans. Signal Proc. 54 (2006): 4634-4643. 

[19] Chen, S. S., D.L. Donoho and M.A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Rev. 43 (2001): 
129-159. 

[20] Cotter, S. E., B.D. Rao, K. Engan and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with 
multiple measurement vectors,” IEEE Trans. Signal Proc. 53 (2005): 2477- 2488. 

[21] Dai, G.-M. and V.N. Mahajan, “Orthonormal polynomials in wavefront analysis: error analysis,” Appl. Opt. 47 
(2008): 3433-3445. 

[22] Dai, W. and O. Milenkovic, “Subspace pursuit for compressive sensing: closing the gap between performance 
and complexity,” IEEE Trans. Inf. Theory55{2009): 2230-2249. 

[23] Davis, G. M., S. Mallat and M. Avellaneda, “Adaptive greedy approximations,” J. Constructive Approx. 13 
(1997), 57-98. 


46 



[24] Daubechies, I. Ten Lectures on Wavelets. (Philadelphia: SIAM, 1992). 

[25] Donoho, D. L., M. Elad and V.N. Temlyakov, “Stable recovery of sparse overcomplete representations in the 
presence of noise,” IEEE Trans. Inform. Theory 52 (2006): 6-18. 

[26] Duarte, M. F. and R.G. Baraniuk, “Spectral compressive sensing,” Appl. Comput. Harmon. Anal. 35 (2013): 
111-129. 

[27] Fannjiang, A., “Compressive imaging of subwavelength structures,” SIAM J. Imag. Sci. 2 (2009): 1277-1291. 

[28] Fannjiang, A., “Compressive inverse scattering I. high-frequency SIMO/MISO and MIMO measurements,” In¬ 
verse Problems 26 (2010): 035008 

[29] Fannjiang, A., “Compressive inverse scattering II. SISO measurements with Born scatterers,” Inverse Problems 
26 (2010): 035009. 

[30] Fannjiang, A., “TV-min and greedy pursuit for constrained joint sparsity and application to inverse scattering,” 
Math. Mech. Complex Syst. 1 (2013): 81-104. 

[31] Fannjiang, A. and W. Liao, “Coherence-pattern guided compressive sensing with unresolved grids,” SIAM J. 
Imaging Sci. 5 (2012): 179-202. 

[32] Fannjiang, A. and W. Liao, “Super-resolution by compressive sensing algorithms,” in IEEE Proc. Asilomar 
conference on signals, systems and computers, 2012 . 

[33] Fannjiang, A., T. Strohmer and P. Yan, “Compressed remote sensing of sparse objects,” SIAM J. Imag. Sci. 3 
(2010): 596-618. 

[34] Fannjiang, C. “Optimal arrays for compressed sensing in snapshot-mode interferometry,” Astron. Astrophys. 559 
(2013): A73-A84. 

[35] Kowalski, M., “Sparse regression using mixed norms,” Appl. Comp. Harm. Anal. 27 (2009): 303-324. 

[36] Mishchenko, M. L, L. D. Travis, and A. A. Lacis. Multiple Scattering of Light by Particles: Radiative Transfer 
and Coherent Backscattering (Cambridge: Cambridge U. Press, 2006). 

[37] Needell, D. and J. A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples”, Appl. 
Comput. Harmon. Anal. 26 (2009): 301-329. 

[38] Pati, Y. C., R. Rezaiifar and P.S. Krishnaprasad, “Orthogonal matching pursuit: recursive function approxima¬ 
tion with applications to wavelet decomposition,” in Proceedings of the 27th Asilomar Conference in Signals, 
Systems and Computers, 1993. 

[39] Rauhut, H. “Stability results for random sampling of sparse trigonometric polynomials,” IEEE Trans. Inform. 
Th. 54 (2008): 5661-5670. 

[40] Rudin, L. and S. Osher, “Total variation based image restoration with free local constraints,” Proc. IEEE ICIP 
1 (1994), 31-35. 

[41] Rudin, L. L, S. Osher and E. Fatemi, ” Nonlinear total variation based noise removal algorithms,” Physica D 
60 (1992): 259-268. 

[42] Tang, C., B. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid”, IEEE Trans- actions on 
Information Theory 59 (2013): 7465-7490. 


47 



[43] Tibshirani, R., “Regression shrinkage and selection via the lasso,” J. Roy. Statist. Soc. Ser. B 58 (1996): 267-288. 

[44] Tropp, J. A., “Greed is good: algorithmic results for sparse approximation,” IEEE Trans. Inform. Theory 50 
(2004): 2231-2242. 

[45] Tropp, J. A., A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. Part I: 
Greedy pursuit,” Signal Process. (Special Issue on Sparse Approximations in Signal and Image Processing) 86 
(2006): 572-588. 

Department of Mathematics, University of California, Davis, CA 95616-8633 
E-mail address: fannjiang@math.ucdavis.edu 


48 



